﻿namespace NCM_MSTest.Alg
{
    /// <summary>
    /// 二元函数插值
    /// </summary>
    public class BinaryFunctionInterpolation
    {
        /// <summary>
        /// 查找最接近u1的u的索引i(后一个点)
        /// </summary>
        /// <param name="n"></param>
        /// <param name="u"></param>
        /// <param name="u1"></param>
        /// <param name="i"></param>
        public static void dij(int n, double[] u, double u1, ref int i)
        {
            int i1 = 1;
            for (i1 = 1; i1 < n; i1++)
            {
                if (u1 <= u[i1])
                    break;
            }
            i = i1;
        }

        public static void LinearIntrp(
            double[] x, double[] y, double[,] f
            , double xi, double yi, ref double fi)
        {
            int n = x.Length;
            int m = y.Length;
            LinearIntrp(n, m, x, y, f, xi, yi, ref fi);
        }

        /// <summary>
        /// 二元分段线性插值
        /// </summary>
        /// <param name="n"></param>
        /// <param name="m"></param>
        /// <param name="x"></param>
        /// <param name="y"></param>
        /// <param name="f"></param>
        /// <param name="xi"></param>
        /// <param name="yi"></param>
        /// <param name="fi"></param>
        public static void LinearIntrp(int n, int m
        , double[] x, double[] y, double[,] f
        , double xi, double yi, ref double fi)
        {
            int i = 0;
            dij(n, x, xi, ref i);
            int j = 0;
            dij(m, y, yi, ref j);

            fi = 0;
            for (int r = i - 1; r <= i; r++)
            {
                for (int s = j - 1; s <= j; s++)
                {
                    double h1 = 1, h2 = 1;
                    for (int k = i - 1; k <= i; k++)
                    {
                        if (k != r)
                            h1 = h1 * (xi - x[k]) / (x[r] - x[k]);
                    }
                    for (int l = j - 1; l <= j; l++)
                    {
                        if (l != s)
                            h2 = h2 * (yi - y[l]) / (y[s] - y[l]);
                    }
                    fi = fi + h1 * h2 * f[r, s];
                }
            }
        }

        /// <summary>
        /// 查找最接近u1的u的索引i（中间点）
        /// </summary>
        /// <param name="n"></param>
        /// <param name="u"></param>
        /// <param name="u1"></param>
        /// <param name="i"></param>
        public static void tij(int n, double[] u, double u1, ref int i)
        {
            i = n - 2;
            for (int i1 = 1; i1 < n; i1++)
            {
                if (u1 <= u[i1])
                {
                    if (i1 != 2
                        && u[i1] - u1 > u1 - u[i1 - 1])
                        i = i1 - 1;
                    else if (i1 == n - 1)
                        i = i1 - 1;
                    else
                        i = i1;
                    break;
                }
            }
        }

        public static void TwiceThreePointIntrp(
            double[] x, double[] y, double[,] f
            , double xi, double yi, ref double fi)
        {
            int n = x.Length;
            int m = y.Length;
            TwiceThreePointIntrp(n, m, x, y, f, xi, yi, ref fi);
        }

        /// <summary>
        /// 二元分段三点插值
        /// </summary>
        /// <param name="n"></param>
        /// <param name="m"></param>
        /// <param name="x"></param>
        /// <param name="y"></param>
        /// <param name="f"></param>
        /// <param name="xi"></param>
        /// <param name="yi"></param>
        /// <param name="fi"></param>
        public static void TwiceThreePointIntrp(int n, int m
            , double[] x, double[] y, double[,] f
            , double xi, double yi, ref double fi)
        {
            int i = 0;
            tij(n, x, xi, ref i);
            int j = 0;
            tij(m, y, yi, ref j);

            System.Diagnostics.Trace.WriteLine($"i-->{i}");
            System.Diagnostics.Trace.WriteLine($"j-->{j}");
            fi = 0;
            for (int r = i - 1; r <= i + 1; r++)
            {
                for (int s = j - 1; s <= j + 1; s++)
                {
                    double h1 = 1, h2 = 1;
                    for (int k = i - 1; k <= i + 1; k++)
                    {
                        if (k != r)
                            h1 = h1 * (xi - x[k]) / (x[r] - x[k]);
                    }
                    for (int l = j - 1; l <= j + 1; l++)
                    {
                        if (l != s)
                            h2 = h2 * (yi - y[l]) / (y[s] - y[l]);
                    }
                    fi = fi + h1 * h2 * f[r, s];
                }
            }
        }

    }
}
